Influence of mass polydispersity on dynamics of simple liquids and colloids 
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We conduct molecular dynamics computer simulations of a system of Lennard- Jones particles, 
polydisperse in both size and mass, at a fixed density and temperature. We test for and quan- 
tify systematic changes in dynamical properties that result from polydispersity, by measuring the 
pair distribution function, diffusion coefficient, velocity autocorrelation function, and non-Gaussian 
parameter, as a function of the degree of polydispersity. Our results elucidate the interpretation 
of experimental studies of collective particle motion in coUoids, and we discuss the implications of 
polydispersity for observations of dynamical heterogeneity, in both simulations of simple liquids and 
colloid experiments. 

PACS numbers: 05.20.Jj, 66.10.-x, 82.70.Dd 



I. INTRODUCTION 

The dynamical behavior of liquids is an area of intense 
current interest. Much of this interest has been moti- 
vated by the desire to understand the progressively slower 
and more complex dynamics of dense, supercooled liquids 
as they are cooled toward the glass transition [|l| . In the 
last few decades, numerous direct insights on dynamical 
motion in liquids have been obtained using molecular dy- 
namics (MD) computer simulations, in which the spatial 
coordinates of particles as a function of time are calcu- 
lated [|| . More recently, experimental studies of colloids 
have used confocal microscopy to track individual parti- 
cles P, 0, pi, pi , thus generating the same type of data on 
microscopic particle motions as is obtained from simula- 
tions. For studying the glass transition, simulations and 
colloid experiments therefore serve as important model 
systems in which the implications of theory can be di- 
rectly tested. 

In both simulations and colloid experiments, fluids 
have been studied in which the particle size is polydis- 
perse. In simulations, size polydispersity is often intro- 
duced to prevent crystallization of the deeply supercooled 
liquid (see e.g. M). In colloid experiments, some degree 
of polydispersity is always present, and depends on the 
process by which colloid particles are produced. To char- 
acterize the polydispersity of colloids, the distribution 
0{(t) of particle diameters a is commonly found (or as- 
sumed) to be Gaussian: 
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where ctq is the average particle diameter, and S charac- 
terizes the width of the distribution M . Polydispersity 
can then be quantified by the value of the dimensionless 
parameter c = S/ctq. 
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In most experimental studies of colloids (see e.g. ||]) a 
system is regarded as efi^ectively monodisperse if c < 0.05. 
For many properties, such as the average liquid structure, 
this is a good approximation. However, dynamical prop- 
erties, especially at the microscopic level, may depend 
sensitively on the nature of microscopic structural fluc- 
tuations, and so may be affected by even small polydis- 
persities. In addition, size polydispersity in real colloids 
leads inevitably to a polydispersity of mass. However, 
some models of polydisperse liquids and colloids consist 
of systems in which particle size varies, but not particle 
mass 1^, 1^ . 

In this paper we seek to isolate and quantify the role of 
size and mass polydispersity on the dynamics of a simple 
liquid system, in particular to assess the need to incor- 
porate mass polydispersity when simulating the dynam- 
ics of realistic systems having size polydispersity. To do 
so, we conduct MD simulations of a system of particles 
interacting via the Lennard- Jones (LJ) potential, poly- 
disperse with respect to both mass and size, as a func- 
tion of c. Our results show that a range of dynamical 
properties (the diffusion coefficient, the velocity autocor- 
relation function, and the non-Gaussian parameter) of 
a polydisperse fluid are systematically shifted from the 
corresponding monodisperse case. We discuss the impli- 
cations of these results for observations of "dynamical 
heterogeneity" in simulations [pi QQ] and in colloid ex- 
periments 



II. POLYDISPERSE LENNARD-JONES LIQUID 

Since our aim is to study generic effects of polydis- 
persity on liquid dynamics, we chose the well-studied LJ 
potential to model interparticle interactions. The LJ po- 
tential is popular for simulations of simple liquids, and 
there exist many studies to which to compare our results. 

In simulations of colloids, the colloidal particles arc of- 
ten modeled as hard spheres, and for many cases this a 
good approximation. However, interactions among col- 
loidal particles can take other forms, and can be explic- 



itly controlled, for example, by attaching soluble poly- 
mer chains by one end to the particle surface to generate 
repulsion, or by adding nonadsorbing soluble polymers 
to the suspension to produce attraction ||ll|. Though 
the present work is motivated by the recent experiments 
studying the dynamics of colloidal particles, we do not 
address the question of how the behavior of a colloidal 
system depends on the shape of the interparticle inter- 
action potential. We also do not take into account the 
influence of a solvent. 

We perform equilibrium molecular dynamics simula- 
tions in three dimensions of a system of A^ = 4000 parti- 
cles interacting via the shifted-force LJ potential, a mod- 
ification of the standard LJ potential. 
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FIG. 1: Effect of polydispersity c on the average liquid struc- 
ture as measured by g{r). 



ij is the potential of interaction of two particles 



Here, Vi 

i and j, separated by a distance r 



e characterizes the 
strength of the pair interaction and is constant for all 
particle pairs. In the shifted-force LJ interaction, the 
LJ potential and force are modified so as to go to zero 
continuously at r = 2.5croj f^nd interactions beyond 2.5cto 
are ignored g . 

Polydispersity is introduced through the particle size: 
CTij ~ {fJi -\- o'j)/2 where cji {aj) characterizes the diam- 
eter of a particle i (j). Particles are assigned cr values 
by random sampling from the Gaussian distribution in 
Eq. nl We also impose a mass polydispersity appropriate 
for the given size polydispersity. The mass of a particle 
i is rrii = mo{ai/ao)^, where ttiq is the mass of a parti- 
cle of size (To- Particle trajectories are evaluated using 
the leap-frog Verlet algorithm Q, using the appropriate 
value of rrii in the equation of motion of each particle. 

Throughout this work we use reduced units. Energy is 
expressed in units of e, length in units of (Tqi the number 
density of particles p in units of ctq" , and temperature T 
in units of s/k, where k is Boltzmann's constant. Time 
t is expressed in units of y/moo^Je. In these units the 
time step used for integrating the particle equations of 
motion is 0.01. 

After equilibration, all quantities are evaluated in the 
microcanonical ensemble. We present data for p = 0.75 
and T = 0.66, a state not far from the triple point of the 
monodisperse LJ fluid (p = 0.85, T = 0.76) ||, |l|. We 
chose this state point so as to avoid the dense, deeply 
supercooled liquid region of the phase diagram of the 
monodisperse LJ system, where spontaneous crystalliza- 
tion could interfere with the evaluation of equilibrium 
properties. We conduct separate simulations for c — 0, 
0.05 and 0.1. 



III. PAIR DISTRIBUTION FUNCTION, 

DIFFUSION COEFFICIENT, AND VELOCITY 

AUTOCORRELATION FUNCTION 



The pair distribution function g{r) that characterizes 
the average liquid structure jl4| is shown in Fig. H for 
each c studied. The effect of increasing polydispersity is 
to reduce the height of, and broaden the peaks associated 
with the successive neighbor shells around each particle. 
However the mean position of each neighbor shell does 
not change noticeably. 

We test for a dependence on c of the bulk transport 
properties by evaluating the diffusion coefhcient D. We 
obtain D from {r^{t)) using the Einstein relation. 
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Fig. H shows the dependence of D on c. We find that at 
fixed p and T, D decreases systematically by about 10% 
as c increases from zero to 0.1. 

Fig. g shows the dependence on c of the velocity auto- 
correlation function ?/j(i) Q: 
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where v(t) is the velocity of a particle at time t. D is 
related to the integral of 4'{'t)i and consistent with the 
decrease of D, the negative part of tlj{t) becomes larger 
in magnitude with increasing c. This trend reflects an 
increase with c of the strength with which single particle 
dynamical properties of the system are retained on a time 
scale comparable to the collision time. 
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FIG. 2: Fractional deviation with c of D relative to Do, its 
value for a perfectly monodisperse system. 



0.04 



0.02 - 



0.00 



-0.02 



-0.04 - 



-0.06 











= 








= 0.05 








= 0.1 


- 










1 /^ ^' ' 

' A/' 


- 








-' 





0.2 0.4 



0.6 
t 



1 1.2 



FIG. 3: Effect of polydispersity c on the velocity autocorre- 
lation function, tp{t). 



IV. NON-GAUSSIAN PARAMETER 

The general non-Gaussian parameter a„(t) is defined 
for integers n > 1 as [na, 



_ (r^"(t)) _ 



(5) 



where c„ = [1 ■ 3 ■ 5 ■ ■ ■ {2n + l)]/3". {r'^"{t)) is the en- 
semble average of the 2n*'' power of the particle displace- 
ments after a time t |lal: 
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Here, ri{t) denotes the position of particle i after a time 
t following a reference time t = in equilibrium. N is 
the total number of particles in the system. 

The functions (r^"(t)) also represent the even moments 
of GsCr, t) , the self-part of the van Hove correlation func- 
tion [Q . For an isotropic fluid made up of particles with 
spherically symmetric interactions, we can restrict our 
attention to Gg {r, t) , the probability density that a par- 
ticle located at the origin at time i = will be found 
within dr of a distance r at time t jlj] : 

Gs{r,t) = {^Y.5{T-\v,{t)^vm\))- (7) 
In terms of Gs{r^t), {r'^^{t)) can be written, 

(r2"(t)) =47r / r^'^Gs{r,t)r^dr. (8) 



In the case of an ideal gas of non-interacting particles 
having a Maxwell-Boltzmann velocity distribution [n4|, 
Gs{r, t) is a Gaussian function of r: 
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where /3 = 1/fcT. In this case, it is readily shown that 
(r2"(i)) = c„(r2(i))" and so a„(i) — 0. For systems 
in which correlations of particle motions are prominent, 
Gs (r, t) is not Gaussian, and the deviation of a„ (i) from 
zero serves to quantify the deviation of Gs {r, t) from the 
Gaussian form. 

In the present work we present results for a2{t), the 
most commonly calculated non-Gaussian parameter: 
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In Fig. ^ we plot a2{t) for three different values of poly- 
dispersity c = 0, 0.05 and 0.1. Qualitatively, there are 
two effects induced by increasing polydispersity: (i) the 
characteristic, intermediate-time peak of 0:2, at approxi- 
mately i = 1, increases in magnitude as c increases; and, 
(ii) the value of a2{t) does not start from zero in the limit 
i — > when c ^ 0. We clarify each of these effects in turn 
below. 

To distinguish the influence of mass and size polydis- 
persity separately, we conduct two new simulations, one 
("size-only") for a system in which the size polydisper- 
sity is c = 0.1, but in which all the particle masses are 
set equal to mg; and another ("mass-only") for which the 
mass polydispersity is taken from our previous "mass- 
and-size" c = 0.1 case, but with all the particle sizes 
then set equal to cto- We compare in Fig. || the resulting 
behavior of a2 as a function of t with the behavior found 
for the monodisperse c = case; and with the case where 
both size and mass are poly disperse with c = 0.1. 

First we focus on the behavior observed near the max- 
imum of a2{t) at approximately t = \. Although the 
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FIG. 4: Variation of a2{t) with c. 



FIG. 6: Log-log plot of 02 as a function of c for a Gaussian 
distribution of particle diameters. 
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FIG. 5: Q2(t) for several types of polydispersity, demonstrat- 
ing that polydispersity of both particle size and mass has a 
greater impact, compared to the monodisperse case, than ei- 
ther size-only or mass-only polydispersity. 



mass-only curve in Fig. || lies above that of the monodis- 
perse system, it still lies well below that corresponding 
to polydispersity of both mass and size. Hence, mass 
polydispersity is not solely responsible for the increase of 
the maximum of oli with c. Interestingly, the size-only 
curve is also well below that corresponding to polydis- 
persity of both mass and size. Even the sum of the devi- 
ations from the monodisperse case of the mass-only and 
size-only curves, is insufficient to account for the height 
of the curve for the system with polydispersity of both 
mass and size. That is, polydispersity of both mass and 
size together has a greater impact on dynamical prop- 
erties at intermediate times than can be obtained from 



polydispersity of mass or size alone. 

Next we turn our attention to the behavior of 02 as 
t — > 0. (For the remainder of this work we will denote 
the limit as i ^- of a2 as a^) In MD simulations of 
a one-component LJ system |l5|, ^ and of a binary LJ 



mixture 



IC 



a2 = 0. However, in the binary LJ mix- 
ture studied, the two species differ in size only and have 
the same mass, in contrast to our system in which the 
masses of particles differ in accordance with the poly- 
dispersity of their sizes. Two of the curves in Fig. 
(c = and "size-only" ) correspond to systems with no 
mass polydispersity, and in both cases a^ = 0. For the 
other two curves (c = 0.1 and "mass-only"), oi°2 adopts 
the same nonzero value. It is clear that the mass poly- 
dispersity is solely responsible for the behavior of o°2 . 

As t —^ 0, the atomic motions in the fluid correspond 
to those of free particles, and the distribution of velocities 
is the Gaussian function given in Eq. |9|. For a monodis- 
perse system, this means that a?, — 0. For a system 
with polydisperse masses, each particle of a given mass 
also samples the Gaussian velocity distribution given in 
Eq. 0. However, the Gaussian distributions sampled will 
have different widths for particles of different m. Conse- 
quently, the form of the total Gg {t, t) function is not in 
general Gaussian because it is a superposition of individ- 
ual Gaussians of different width. The result is a non-zero 
value of a2, as found in our simulations. 

Since the limit i ^ corresponds to the free-particle 
limit for atomic motion, we can calculate the non- 
Gaussian parameter for a polydisperse system as i ^ 
exactly. Consider a system of N particles in which there 
are M species (labeled by index j) each having Nj par- 
ticles of mass ruj. The moments (r^"(t)) can be found 
using a modified form of Eq. H appropriate for an M- 



component system: 
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Eq. |ll| can be rewritten as, 
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where fj — Nj/N is the fraction of particles of species j 
and 
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where the sum is over particles only of species j. For each 
species, the atomic motion as t ^ is also described by 
Eq. P with the appropriate value of ttt, = toj, and hence, 
the moments {r'^{t))j and {r^{t))j can be found in the 
limit i — > by substituting Eq. |9 into Eq. for each j. 
The result is. 
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constant for all choices of p, T and interparticle interac- 
tion. 

We apply the above result to the case of the Gaus- 
sian distribution of particle diameters given in Eq. ^. We 
substitute 9 = Oq with ctq = 1 in Eq. |l^ and calculate 
a2 as a function of c (Fig. O). We evaluate the integrals 
in Eq. |l8| numerically, replacing the limits of integration 
(0,oo) with [0. 01(70, 2cro]. This avoids the divergence of 
the integrands at cr = 0, and in any case is more physi- 
cal, since a real distribution of particle sizes would have 
a non-zero lower bound, and a finite upper bound. As 
seen in Fig. 0, aj oc c^ for small c, but increases more 
rapidly than this for c > 0.1 The predictions of Eq. |6| are 
in agreement with our simulation results. For c = 0.05 
Eq. na gives a2 = 0.02344, while our simulation gives 



a^ = 0.027; 



for c = 0.1 we obtain a^ = 0.10781 and 



0.113, respectively. 

We can also use Eq. |lj to calculate ctj for a system 
in which the distribution of sizes is not Gaussian. For 
example. Sear [Q simulated a system of hard spheres, 
polydisperse in both size and mass, using a "hat" func- 
tion of width w: 0(a) = {wao)~^ for cro(l ~ w/2) < cr < 
cro(l + u'/2), and 9{(t) — otherwise, and m ^ a^ . For 
this case, we are able to solve Eq. 08 exactly, giving. 
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(19) 



For w = 0.3, Eq. |l| gives a.2 = 0.06916, while for uj = 0.7, 
cti — QA2T1, in agreement with the simulation results in 
Fig. 6of Ref. pi. 



The value of a^ for the multicomponent system can 
then be found by using Eqs. n3, nJ and na in Eq. nfl: 
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This expression highlights that a2 may not equal zero for 
a system with polydisperse masses. 

If the polydispersity is expressed as a continuous dis- 
tribution of masses 4>{m), Eq. O generalizes to, 
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Experimental studies of colloids typically characterize 
polydispersity not in terms of the masses, but in terms 
of the particle diameters, described by 0(a). Assuming 
that the particle mass m is proportional to a^, and that 
0(to) dm — 9{a) da, Eq. H^ becomes, 



J,^a-Hia)da 
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Note that the value of a2 therefore depends only on the 
shape of the mass distribution function, and is otherwise 



V. DISCUSSION 

When a2 has been extracted via confocal microscopy 
in colloid experiments, values as high as a2 ~ 1.5 have 
been observed |^, |j, g . In these studies the polydispersity 
ranged from c = 0.01 to 0.1. Kasper, et al. observed 
that a2 is not zero for all values of the volume fraction oc- 
cupied by the colloidal particles. Marcus, et al. [W found 
that a2 is nonzero and increases with increasing volume 
fraction. Weeks, et al. [^ found that of aj is approxi- 
mately constant for small volume fractions but grows for 
higher volume fractions. In general, a2 was found to in- 
crease with the volume fraction occupied by the colloid 
particles, in contrast to the absence of any density depen- 
dence in Eq. 08. In the case of real colloids, the behavior 
of a2 as f — > is complicated by the fact that solvent- 
induced hydrodynamic forces among particles potentially 
introduce strong, short-time-scale correlations of particle 
velocities, invalidating the free-particle assumption that 
is the basis of Eq. ||. The large difference between the be- 
havior of Qfj found for these systems, and that predicted 
by Eq. |l^ demonstrates that polydispersity alone cannot 
account for the observed values of aj , and that hydrody- 
namic effects indeed dominate the short-time dynamical 
behavior of real colloids. 



At intermediate times, we find that the peak value of 
a2 increases with c. The maximum value of a2 has been 
shown [Q to correlate to the degree of dynamical het- 
erogeneity present in the system: that is, transient, spa- 
tially correlated groups of particles whose characteristic 
structural relaxation time differs from the mean. Our 
results therefore suggest that dynamical heterogeneity, 
prominently observed at lower T and higher p than stud- 
ied here, may be enhanced as polydispersity increases. 
One source of this enhancement may be the influence of 
mass polydispersity on spontaneously occurring density 
fluctuations, that in turn control the development of dy- 
namical heterogeneities. In general, the occurrence and 
quantification of dynamical heterogeneity in a polydis- 
perse system is likely to be more complicated than in a 
monodisperse (or even bi-disperse) system. At the same 
time, since we observe a slowing of the dynamics with in- 
creasing polydispersity, the dynamical heterogeneity may 
be more prominent and longer-lived in a polydisperse sys- 
tem, and so may facilitate the study of these complex 
structures. 

In summary, we have illustrated that a system of par- 
ticles with polydispersity of both mass and size is a more 
realistic model (compared to models without mass poly- 
dispersity) for studying the dynamics of colloids in MD 
simulations. Our results show that typical polydisper- 



sities found in real systems can induce an influence of 
the order of 10% on dynamical properties. This is a 
small effect for studies, such as those near a glass tran- 
sition, where relaxation times may vary by several or- 
ders of magnitude. At the same time, knowledge of the 
amount and direction of the impact of polydispersity on 
dynamics is required because polydispersity is so com- 
monly found in systems studied both in simulations and 
experiments. This knowledge is also crucial for precise 
tests of theories, particularly those formulated for per- 
fectly monodisperse systems. We also note that our re- 
sults can be tested experimentally in colloids by deliber- 
ately varying the polydispersity of the studied colloidal 
system. Indeed, by varying c alone it might be possible 
to study an "isothermal-isochoric glass transition" (i.e. a 
glass transition at both fixed T and p) by setting up an 
appropriate series of colloidal systems where the polydis- 
persity is progressively increased. 
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